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Abstract 



. We describe a method for restoring information lost during statistical thinning in extensive air shower simulations. By converting 
Vveighted particles from thinned simulations to swarms of particles with similar characteristics, we obtain a result that is essentially 



identical to the thinned shower, and which is very similar to non-thinned simulations of showers. We call this method dethinning. 
Using non-thinned showers on a large scale is impossible because of unrealistic CPU time requirements, but with thinned showers 
that have been dethinned, it is possible to carry out large-scale simulation studies of the detector response for ultra-high energy 
^) cosmic ray surface arrays. The dethinning method is described in detail and comparisons are presented with parent thinned showers 
and with non-thinned showers. 
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^ 1. Introduction 
HH ! 

■ In the study of ultrahigh energy cosmic rays (UHECR) two 
Q-imain experimental techniques have been used: detection of the 
Q fluorescence light emitted by nitrogen molecules excited by the 
^ passage of particles in extensive air showers, and detection of 
c/2 the particles themselves when they strike the ground by deploy- 
I ^ I ing an array of particle detectors over a large area. Analysis of 
data from a fluorescence detector involves making a detailed 
Monte Carlo simulation of the shower, atmosphere, and detec- 
tor Ol 01 ■ Only by this technique can the aperture be calcu- 
lated as a function of energy. Simulation of a large number of 
complete showers can not be performed using programs like 
Cn CORSIKA d or AIRES (3] because the CPU requirements are 
too large. The approximation technique called thinning is there- 
fore used, in which particles are removed from consideration in 
the shower generation and other particles in similar regions of 
phase space are given weights to account for the loss. Since a 
fluorescence detector is sensitive to the core region of a shower, 
.where 10'^ charged particles occur at shower maximum for a 
■10^" eV event, thinning does not harm the accuracy of the sim- 
ulation. However, for an array of surface detectors (SD), where 
' several km from the core the density of particles is low, the thin- 
ning approximation produces an unacceptably coarse simula- 
tion of a shower The average density of particles, as a function 
of radius from the core, in a thinned shower is approximately 
correct, but the fluctuations about the average are much larger 
than in a shower generated without using the thinning approxi- 
mation. The result of this is that the Monte Carlo technique has 
been available to those analyzing SD data only in a very limited 
way. 
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We have developed a method of performing an accurate 
Monte Carlo simulation of the surface detector of the Telescope 
Array (TA) experiment. This method consists of three parts: (1) 
generating 100 non-thinned CORSIKA showers above 10'^'^ 
eV using many computer nodes operating in parallel jsl, (2) 
using these showers to develop a method (called dethinning) of 
replacing much of the shower information lost in thinning, and 
(3) generating a large number of dethinned showers, including 
a detailed simulation of the TA surface detector performance, 
and comparing histograms of important quantities between the 
data and the Monte Carlo simulation. Reference |5] describes 
a method for generating CORSIKA showers using many com- 
puting nodes in parallel. The present work describes the second 
step of the method: dethinning. A future paper will describe 
the third step. This method has proved quite successful, and 
has allowed us to calculate the aperture of the TA surface de- 
tector even in the energy range where its trigger efficiency is not 
100%. 

The idea of replacing the information lost in thinning was 
first introduced by P. Billoir [6]. The basic idea in reference 
f^"] and in our work is the same: start with a thinned shower, 
maintain the average density of particles, and smooth the distri- 
bution to get the correct amount of fluctuations. He considered 
CORSIKA output particles striking the ground in the vicinity 
of a surface detector (specifically a detector of the Pierre Auger 
experiment), and, by an oversampling technique based on the 
CORSIKA output particles, predicted what that detector should 
observe. He named this technique unthinning. He continued 
on to estimate several systematic biases to which his method 
might be sensitive, by studying the properties of thinned and 
unthinned showers. 

Our method, as described in this paper, differs in two ways. 
First, as the dethinning process we take each CORSIKA out- 
put particle with weight w and from it generate a swarm of w 



Preprint submitted to Astroparticle Physics 



January 19, 2013 



particles. We perform this step for the entire set of CORSIKA 
output particles. The details of this generation matter, and are 
described in this paper Second, we compare the resulting de- 
thinned showers with showers of almost identical characteris- 
tics generated without using the thinning approximation. This 
is a direct way of testing the accuracy of the dethinning pro- 
cess. Since real surface detectors for ultrahigh energy cosmic 
rays sample showers coarsely, and measure only the time distri- 
bution of the number of particles that strike them, this is the im- 
portant aspect of the comparison process. We present our com- 
parisons between dethinned showers and non-thinned showers 
in this manner. 

2. Dethinning Description 

In a thinned EAS simulation, particles are discarded from the 
simulation in order to conserve computation time. In the case 
of CORSIKA for a given thinning level, e,^, if the energy 
sum of all j secondary particles falls below the thinning energy 



£,,,£■0 > y^jEj. 



(1) 



then only one randomly assigned secondary particle survives 
with probability 



(2) 



If the energy sum is greater than the thinning energy, then sec- 
ondary particles with energy below the thinning energy survive 
with probability 



Pi = Ei/(s,i,Eo). 



(3) 



In both cases, surviving particles have their weight multiplied 
by a factor of w, = 1 / pj. Thus the weight of a particle reaching 
the end of the simulation after passing through k vertices is 



(4) 



For sufficiently low values of e,/,, it is clear that the thinned 
simulation output can be thought of as an accurate sample of 
secondary particle types, trajectories, and positions compared 
to a non-thinned simulation, for the observable parts of the 
shower. In this situation, for a particle of weight, w,-, the sim- 
ulation, on average, removed w, - 1 particles from a similar 
position in phase space. (Of course, if the value of e,/, is in- 
creased, this situation is no longer valid because the thinned 
simulation no longer has the same distribution of particle types, 
trajectories, and energies as the parent shower) By comparing 
a dethinned shower with a similar non-thinned shower one can 
determine the accuracy of the sampling. The pivotal questions 
are then: (1) How can a thinned sample be used to reconstruct 
the fuU simulation? and (2) What is the maximum value of s,i, 
for which the thinned sample accurately represents the parent 
shower's particle types, etc.? 

We address the first question by describing the process by 
which we dethin the showers. The original CORSIKA shower 
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Figure 1: Geometry for a "Gaussian cone" with a vertex placed at arbitrary 
position on tlie trajectory of the weighted particle 



consists of a list of output particles (plus their weights, types, 
energies, positions, angles, and arrival time) that have struck 
the ground. Dethinning consists of adding particles to this list. 
For every CORSIKA output particle of weight w we add w - 1 
particles to the list. When this is completed the weight of each 
particle is set to 1 . To insert these particles we use the following 
procedure (see Figure [T}. 

1. Choose a vertex point on the trajectory of the weighted 
particle, in the way given in the next paragraph. 

2. Choose a point in a cone centered on the output particle's 
trajectory, weighted by a two-dimensional Gaussian dis- 
tribution with a sigma of a few degrees (as described in 
Section 3). This will be the inserted particle's trajectory. 

3 . Project the inserted particle to ground level, assign it a time 
and energy (as described in Section 3), and add it to the 
particle list of the dethinned shower 

4. Perform steps 2-3 w-l times. For the case where w is not 
an integer, add one particle randomly based on the decimal 
part of w. 

There is a maximum distance from the ground that one can 
choose for the vertex in item 1 above, which is set by the re- 
quirement that no particle can have an arrival time that precedes 
the arrival of the shower front. A too-early arrival time occurs 
when the total time-of-flight from the point of first interaction, 
xo, to the vertex point and then to the position on the ground 
of the generated particle is less than the time-of-flight directly 
from Xo to final particle position. This can be corrected by fix- 
ing the position of the vertex point to a position where the time- 
of-flight from Xo to the imaginary vertex and then onward to the 
final position of the weighted particle, x, is equal to the diff'er- 
ence in the arrival time of the weighted particle, f,, and the time 
of first interaction, to. This condition is: the distance along the 
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Figure 2: In order to ensure temporal consistency in tlie EAS simulation, we 
require tj - to > d/j + df^ , where tj is the recorded anival time for weighted 
particle and Iq is the time of first interaction. 

weighted particle trajectory, pi, between Xi and the vertex point 
is 

_ c\ti - tof - |Xi - XqI^ 
2(c(ti - to) - (Xi - Xo) ■ Pi) 

where c is the speed of Hght. We should emphasize that D,„ax 
is the maximum separation between the vertex point and the 
ground. Any shorter separation will also generate dethinned 
particles that are temporally consistent. The calculation of D„,ax 
is graphically illustrated in Figure |2] 

The second important question pertains to selecting a value 
of Srh- It should be recalled that in order for this algorithm to 
work, Sfh must be sufficiently small so that the thinned simu- 
lation contains a large enough sample of output particles that 
the distribution of particle types, trajectories, and positions is 
the same as in the non-thinned simulation. We addressed this 
question phenomenologically by comparing non-thinned and 
dethinned showers. Our conclusion was that for e,/, = 10"^, this 
algorithm could be applied without further adjustment. Con- 
versely, for e,;, - 10"^, the thinned sample did not contain 
enough output particles to successfully mimic the properties 
of a non-thinned shower. The option in between, e,/, = 10"^, 
proved to be a borderline case. We found that dethinning could 
be successfully implemented by careful adjustment of the avail- 
able parameters. Because 10"* showers take 1/10 the time to 
generate as 10"^ showers, and require 1/10 as much storage 
space, Sth - 10"* proved to be highly desirable so we chose 
to focus our efforts there. 

3. Adjusting the Dethinning Algorithm 

In adjusting the dethinning algorithm, we have sought agree- 
ment between dethinned and non-thinned simulations for all 
measures relevant to observation by the TA surface array. 
These measures include: distributions of secondary particle 
position and time, particle type, incident angle, and energy. 



These measures were tested by comparing thinned versus de- 
thinned versions of the same simulation and then subsequently 
comparing dethinned simulations versus non-thinned simula- 
tions with identical input parameters. Tuning was not consid- 
ered complete until all measures agreed for lateral distances 
[500, 4500] m from the shower core. 

In the first step, secondary particle spectra for thinned and 
dethinned versions of the same shower were compared with re- 
spect to particle type (photons, electrons, and muons), incident 
angle with respect to the ground, and position within the shower 
footprint. The algorithm is tuned so that the particle fluxes af- 
ter dethinning were consistent with those seen in the original 
thinned output. 

In the second step, distributions of particle fluxes over 
6 X 6 m^ detector-size areas are compared for dethinned and 
non-thinned simulations. For this purpose, a library of more 
than 100 non-thinned showers was generated with parallelized 
CORSIKA 1 5]. This library contains showers with Eq = 
[10"*■^ 10''^-^]eV, Oo = [0,60]°, and proton and iron primary 
particle types. When identical input parameters are used for 
thinned and non-thinned simulations, the resulting simulations 
are not identical. It is therefore necessary to normalize the net 
secondary particle fluxes of the non-thinned simulation with 
respect to the thinned simulation. Once this normalization is 
accomplished the dethinning algorithm can be further refined 
so that shower particle fluctuations are consistent between de- 
thinned and non-thinned simulations. A further check is done 
by simultaneously examining dethinned versus non-thinned 
comparisons over many simulations without normalizing the 
non-thinned showers. By utilizing the 100 non-thinned showers 
for this comparison, we ensure that the thinning-dethinning pro- 
cess does not bias the energy scale with respect to non-thinned 
showers. 

The result of adjusting the parameters of the dethinning pro- 
cess is as follows: 

1 . Angle subtended by Gaussian cone: Set to /3d where d is 
the lateral distance from the shower core for the weighted 
particle and /3 = 3°/km for electromagnetic particles and 
l°/km for muons and hadrons. The values of yS are the 
minimum necessary to dethin simulations with e,/, = 10"*. 
A smaller value of e,/, enables the use of smaller /? values. 

2. Energy perturbation: Vary the energy of each particle in 
swarm about a ±10% fractional Gaussian distribution cen- 
tered on the energy of the original particle. This correction 
smooths the secondary particle spectra. 

3. Minimum lateral distance: Since any detector within a few 
hundred meters of the shower core is saturated by the large 
number of particles in that part of a shower, it is not nec- 
essary to simulate the central part of a shower; not doing 
so saves a great deal of CPU time as well. We therefore 
do not consider in the dethinning process any CORSIKA 
output particle within a distance r„„„ of the core, and to 
avoid biases retain only particles farther than r'^^^^^^ from 
the core. For the case of - 10"*, r„„„ > 100 m, and 

• - f,nin ^ 200 m. 
mm 

4. Particle acceptance: Some particles in the swarm with 
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longer trajectories than the original weighted particle, if 
they were followed in the simulation, would not reach the 
ground. We therefore introduce an acceptance for parti- 
cles in the swarm with probability: P - e^^'^^'^, where Ax 
is the difference in slant depth between the trajectories and 
e = 50 g/cm^. This is an appreciable correction only for 
showers with large zenith angle. 
5. Height of Gaussian cone: When one sets the vertex a 
height AiM.v above the ground, this aligns the time at the 
ground for reintroduced particles, but muons that are cre- 
ated late in the shower then have too wide a spatial distri- 
bution A solution is to set the vertex distance to the smaller 
of Dmax and 

D' = |xi-xo|-X-'(xi,xo,a/!), (6) 

where xq is the point of first interaction, h is the gener- 
ation of the hadron from which the particle originated, 
a = 30 g/cm^, and X"'(xi,xo, a/;) is the distance equiv- 
alent of ah slant depth on the trajectory from xo to Xj. 

4. Comparing Simulations: An Example 

In Section [3] the method used to tune the dethinning algo- 
rithm for E,ii - 10 thinned showers was described. We now 
consider two examples of comparisons between the dethinned 
result and the parent thinned shower and with a very simi- 
lar non-thinned shower For the comparisons, we use COR- 
SIKA v6.960 High energy hadronic interactions are mod- 
eled by QGSJET-II-03 10], low energy hadronic interactions are 
modeled by FLUKA2008.3c f^|?], and electromagnetic inter- 
actions are modeled by EGS4 [ 101 . For the thinned shower, 
s - 10"^ and additionally, we apply the thinning optimization 
scheme proposed by Kobal (vi\. 

For both comparisons, the shower footprint is divided into 
eight 500 m thick ring-shaped segments from 500 to 4500 m 
in lateral distance. Each lateral ring is further divided into six 
pie-shaped wedges with respect to the rotation angle about the 
shower axis. 

For the first comparison, we divide the particle flux into ten 
cos di =0.1 bins, where 0, is the incident angle of an individ- 
ual particle with respect to the ground. Three particle types are 
considered: photons, electrons, and muons (all other types are 
relatively scarce). Each bin is then histogrammed with respect 
to energy. This results in 8x6x 10x3 - 1440 discrete secondary 
particle energy spectra. By scanning through these spectra, any 
discrepancies in particle flux generated by dethinning can be 
readily identified. Figures |3] and |4] show examples, typical of 
all 1440 spectra, of these comparisons. The agreement is excel- 
lent, showing that the process of smoothing the distribution of 
particles does not change the angular or energy distribution of 
particles from the parent thinned shower 

Having established that dethinning maintains the large-scale 
secondary particle fluxes from parent thinned simulations, we 
then turn to comparisons with non-thinned showers produced 
with a parallelized version of CORSIKA |5]. Because it 
is structurally impossible in CORSIKA to produce identical 
thinned and non-thinned simulations, for this comparison it is 



necessary to normalize the net particle flux of the non-thinned 
and thinned (not dethinned) simulations. This is done sepa- 
rately for each wedge-shaped region of the shower footprint and 
each particle type. 

For the comparison, we consider the same segments in the 
shower footprint as for Figures |3]and 2] The segment is then 
further divided into 6 x 6 m^ tiles covering the distances from 
500 m to 4500 m from the core. These tiles are then projected 
onto the ground, and for each tile we tabulate the time, fi/io, 
when 10% of the total particle flux has arrived, the time, ti/2, 
when 50% of the total particle flux has arrived, and the flux of 
all photons, electrons, and muons. The times are then corrected 
for the time offset between the positions of each segment on the 
ground and in the plane normal to the EAS. Figures |5] through 
|9] show the results of this comparison. They show that for the 
dethinned shower the time of arrival and number of particles 
per tile agree very well with that of the non-thinned shower. 

5. Conclusion 

The aim of this dethinning method is to use a thinned simula- 
tion to reconstruct, on a statistical basis, information lost in the 
thinning process. Since thinning preserves mean particle den- 
sities as a function of radius from the shower core, but intro- 
duces large biases into the distribution of the densities (e.g., the 
RMS of the particle distribution), dethinning is designed to be a 
smoothing procedure. The dethinning process is tuned to main- 
tain mean particle densities from parent thinned showers, and 
reproduces distributions of arrival times and numbers of parti- 
cles striking counter-size areas in non-thinned showers. Since 
these are the distributions to which surface detectors of experi- 
ments are sensitive, dethinned showers can be used in place of 
non-thinned showers in comparisons with data. 

This method has three primary limitations. We require that 
s,h < 10"^ and lateral distances be restricted to less than 4500 m 
from the shower core. Beyond these limits we cannot reliably 
control for artificial fluctuations. We have tested this process for 
9q < 60° but have not yet examined the case of more inclined 
showers. 

We have compared dethinned shower simulations against 
both thinned and non-thinned simulations, and shown that 
the dethinning process reproduces the characteristics of 
CORSIKA/QGSJET-II-03 showers generated wifliout thinning. 
In a future paper in this series we will show that the resulting 
showers reproduce the characteristics of the TA surface array 
data. 

Dethinning is proving to be a powerful tool for studying and 
understanding UHECR observations by enabling a thorough 
simulation of surface array data. This enables a direct com- 
parison between Monte Carlo simulations and TA surface array 
data and results in a more complete understanding of the re- 
sponse of the detector This understanding leads to the ability to 
accurately assess detector aperture for efficiencies far less than 
100%, which in turn leads to significant improvements in the 
measurement of the cosmic ray spectrum, and in the estimation 
of the detector exposure to the sky for astronomical studies . 
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Figure 3: Comparison of secondary particle spectra with and without dethinning for a thinned simulation of a protonic EAS with primary energy Eg = lO''"' eV and 
primary zenith angle 60 = 45°. In both cases, the secondary particles whose ground position was within a region enclosed by shower rotation angles, <!> = [-30°, 30°] 
(with respect to the primary azrmuthal direction) and lateral distances, r = [500m, 1000m] were tabulated according to particle type, incident angle with respect to 
the ground, 0,-, and kinetic energy. In the case of the thinned simulation, each secondary particle with weight, w,-, was treated as w; identical particles. The resulting 
spectral comparisons are shown in cos 0, = 0.1 increment bins for a) photons, b) electrons, and c) muons. For each histogram, good agreement is observed between 
thinned simulations (gray) and dethinned (black). 
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Figure 4: Comparison of secondary particle spectra with and without dethinning for a thinned simulation of a protonic EAS with primary energy Eg = 10 eV 
and primary zenith angle do = 45°. In both cases, the secondary particles whose ground position was within a region enclosed by shower rotation angles, 
<J> = [150°, 210°] (with respect to the primary azimuthal direction) and lateral distances, r = [1500m, 2000m] were tabulated according to particle type, incident 
angle with respect to the ground, 9,, and kinetic energy. In the case of the thinned simulation, each secondary particle with weight, w;, was treated as w; identical 
particles. The resulting spectral comparisons are shown in cos 6; = 0.1 increment bins for a) photons, b) electrons, and c) muons. For each histogram, good 
agreement is observed between thiimed simulations (gray) and dethinned (black). 
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Figure 5: Comparison of the distribution of rise times, ?i/io, wliere 10% of tlie total particle flux has arrived for a given 6 X 6 m^ segment in plane normal to shower 
trajectory for lO"" eV protonic EAS simulations with a primary zenith angle of 45°. For this comparison, ?i/io was measured for segments within a) a region 
enclosed by shower rotation angles, (J = [-30°, 30°] (with respect to the primary azimuthal direction) and lateral distances, r = [500m, 1000m] and b) a region 
enclosed by shower rotation angles, <b = [150°, 210°] (with respect to the primary azimuthal direction) and lateral distances, r = [1500m, 2000m]. In both cases, 
the distribution of /|/io values is consistent for the dethinned (black) and non-thinned (blue) simulations while the thinned (red) simulation is quite different. 
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Figure 6: Comparison of the distribution of median arrival times, ?i/2, where 50% of the total particle flux has arrived for a given 6 X 6 m^ segment in plane normal 
to shower trajectory for lO"" eV protonic EAS simulations with a primary zenith angle of 45°. For this comparison, ti/i was measured for segments within a) a 
region enclosed by shower rotation angles, (D = [-30°, 30°] (with respect to the primary azimuthal direction) and lateral distances, r = [500m, lOOOm] and b) a 
region enclosed by shower rotation angles, <D = [150°, 210°] (with respect to the primary azimuthal direction) and lateral distances, r = [1500m, 2000m]. In both 
cases, the distribution of fi/2 values is consistent for the dethinned (black) and non-thinned (blue) simulations while the thinned (red) simulation is quite different. 
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Figure 7: Comparison of the distribution of plioton flux measurements for 6 X 6 m- segments in plane normal to shower trajectory for lO"" eV protonic EAS 
simulations with a primary zenith angle of 45°. For this compaiison, photon flux measurements were done for segments within a) a region enclosed by shower 
rotation angles, d) = [—30°, 30°] (with respect to the primary azimuthal direction) and lateral distances, r = [500m, 1000m] and b) a region enclosed by shower 
rotation angles, <D = [150°, 210°] (with respect to the primary azimuthal direction) and lateral distances, r = [1500m, 2000m]. In both cases, the distribution of 
photon flux values is consistent for the dethinned (black) and non-thinned (blue) simulations while the thinned (red) simulation is quite different. 
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Figure 8: Comparison of the distribution of electron flux measurements for 6 X 6 m^ segments in plane normal to shower trajectory for lO"" eV protonic EAS 
simulations with a primary zenith angle of 45°. For this comparison, electron flux measurements were done for segments within a) a region enclosed by shower 
rotation angles, d) = [—30°, 30°] (with respect to the primary azimuthal direction) and lateral distances, r = [500m, 1000m] and b) a region enclosed by shower 
rotation angles, <D = [150°, 210°] (with respect to the primary azimuthal direction) and lateral distances, r = [1500m, 2000m]. In both cases, the distribution of 
electron flux values is consistent for the dethinned (black) and non-thinned (blue) simulations while the thinned (red) simulation is quite difi'erent. 
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Figure 9: Comparison of the distribution of muon flux measurements for 6 X 6 m^ segments in plane normal to shower trajectory for lO"" eV protonic EAS 
simulations with a primary zenith angle of 45°. For this comparison, muon flux measurements were done for segments within a) a region enclosed by shower 
rotation angles, (D = [-30°, 30°] (with respect to the primary azimuthal direction) and lateral distances, r = [500m, 1000m] and b) a region enclosed by shower 
rotation angles, <D = [150°, 210°] (with respect to the primary azimuthal direction) and lateral distances, r = [1500m, 2000m]. In both cases, the distribution of 
muon flux values is consistent for the dethinned (black) and non-thinned (blue) simulations while the thinned (red) simulation is quite dift'erent. 
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